Real-Time-RG Analysis of the Dynamics of the Spin-Boson Model 
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Using a real-time renormalization group method we determine the complete dynamics of the spin- 
boson model with ohmic dissipation for coupling strengths a < 0.1 — 0.2. We calculate the relaxation 
and dephasing time, the static susceptibility and correlation functions. Our results are consistent 
with quantum Monte Carlo simulations and the Shiba relation. We present for the first time reliable 
results for finite cutofi^ and finite bias in a regime where perturbation theory in a or in tunneling 
breaks down. Furthermore, an unambigious comparism to results from the Kondo model is achieved. 
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Introduction. The spin-boson model is one of the most 
fundamental quantum dissipative systems It plays 

an important role in describing defect tunneling in solids 
|B,Q, quantum tunneling between flux states in a SQUID 
|5|, electron tunneling between quantum dots or elec- 
tron transfer in chemical and biological reactions . The 
model provides a nontrivial description of dissipation in a 
quantum system and has attracted much interest due to 
its simplicity and applicability. It consists of a two-state 
(spin) system coupled linearly to an infinite-dimensional 
harmonic oscillator (boson) bath with the Hamiltonian 
H = Ha + Hb + V, where Ha = -fcr^^ + f cr^, Hb = 
J2q ^^qa\aq, and V = ^-Y,q 9q{a\ + aq). Here, and 
are the usual Pauli matrices. A is the tunnel matrix ele- 
ment and e the energy difference between the two states. 



model on the anisotropic Kondo model ||T^, and solve 
the latter exactly using numerical renormalization group 
(NRG) 1 14], Bethe ansatz ||l^, or conformal field theory 
(CFT) [16|. However, NRG and Bethe ansatz provide 
only spectral properties or dynamics at very short time 
scales [0, and CFT has solved so far only the unbiased 
case e = for the diagonal elements of p{t). Further- 
more, and most importantly, the mapping on the Kondo 
model can not be proven rigorously, and the relation of 
the parameters is not precisely known P,p^ . It is known 
that the mapping is incorrect for finite cutoff D, but it 
is at least established that the scaling behaviour agrees 
with that of the spin- boson model [|l|,0[6). 

In this paper, we will present for the first time a solu- 
tion of the complete dynamics of the spin-boson model 
(flg) creates (annihilates) a boson with energy huq. for a ^ 0.1 — 0.2. We study the diagonal and nondi- 

agonal part of the reduced density matrix p(t) starting 
from an arbitrary nonequilibrium state pa. From the 
asymptotic behaviour we determine the relaxation and 
dephasing time. Furthermore, we calculate the spin sus- 
ceptibility and correlation functions. Especially, and in 
contrast to many other methods, we solve directly the 
spin-boson model and present results at finite cutoff D 
and finite bias e. Therefore, our results provide for the 
first time the possibility for a quantitative and unambi- 
gious comparism to results obtained from the mapping 
of the spin-boson model on the Kondo model. We find 
that for e = the relaxation parameters agree with those 
of CFT and in the scaling limit the susceptibility agrees 
rather well with Bethe-ansatz results, but for the lat- 
ter deviations occur at finite bias. To demonstrate the 
reliability of our results, we show the consistency with 
chromostochastic quantum dynamics (CSQD) ||l^, and 
check the Shiba-relation as well as the scaling behaviour. 

To obtain our results, we will use a recently developed 
real-time renormalization group (RTRG) method ||l9|,^ . 
This method has been successfully applied to equilibrium 
problems ]l9| , pl[, a nd to the study of nonequilibrium sta- 
tionary states |2^. Here we will apply it for the first 



The quantity gq represents the coupling strength of the 
two-state system to the coordinate of the q*'* oscillator. 
The coupling to the environment is completely defined 
by the spectral density J{uj) — T^J2q9qH'^ ~ ^q)-< which 
is usually parametrized by J{uj) = 27raa;"+^e^'^/^. The 
case n — corresponds to the Ohmic bath, which we 
want to consider in the following. Therefore the coupling 
to the bath is now characterized by the coupling strength 
a and the high-energy cutoff D. 

Except for some special parameter values an exact so- 
lution of the spin-boson model is not known ^ . Most of 
the studies are based on perturbative approaches in a j^] 
or A 1^. The latter is known as the noninteracting blip 
approximation (NIBA) and gives very reliable results at 
zero bias for the diagonal matrix elements of the reduced 
density matrix p{t) of the two- level system. Recent real- 
time quantum Monte Carlo (QMC) simulations [ p^JTT| 
provide also reliable information on the nondiagonal ele- 
ments of p{t) , but they were only applied to the symmet- 
ric case and the correct long-time behaviour has not been 
checked. Flow equation methods based on infinitesimal 
unitary transformations JT2| have reproduced the Shiba- 
relation up to a ^ 0.025 — 0.05 (with an error of 3 — 10%) , 
but they only addressed spectral properties of the system. 
Other approaches try to use a mapping of the spin-boson 



time to the dynamics of the reduced density matrix, and 
generalize it to the evaluation of correlation functions as 
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well. Since this technique is rather straightforward and 
easily generalized to other models, this may open a new 
possibility for the study of various kinds of dissipative 
quantum systems, like e.g. many-level systems, magnetic 
nanoparticles interacting with phonons, coupled quan- 
tum dots or other kinds of dissipative environments. 

Kinetic equation and RTRG approach. The RTRG ap- 
proach is based on a kinetic equation for the density ma- 
trix, see Ref. for details. We only mention the main 
steps: 

1 . The time evolution of the reduced density matrix is 
written in Liouville space as p{t) = Tr^ exp {—iLt)poPg, 
where Tr^ denotes the trace over the bath degrees of 
freedom, and L = [H, ■] = Lq + Lb + Ly is the Liouville 
operator. The initial density matrix is assumed to de- 
couple into an arbitrary nonequilibrium distribution 
for the two-level system and an equilibrium distribution 
p^g for the oscillator bath. 

2. The propagator exp (—iLt) is expanded in the in- 
teraction part Lv, and the trace Tr^ over the bath de- 
grees of freedom is performed by application of Wick's 
theorem. In this way one obtains a series of terms where 
vertices C of the local system (originating from the (7^/2 
factor in V) are connected by pair contractions (t) of 
the bath. Here, p = ± indicates wether the interaction 
takes place on the forward or the backward propagator. 
With j = J2q9q{al + aq) and 7(t) = TrBPBj{t)j, we 
obtain -/Pp' (t) = R{t) + ip'S(t) with -f{t) = R{t) + iS{t) 
and 



R{t) = -2aRe{{TrTf 



S{t) = -2alm{ 



sinli (7rr(t-i/D)) 



{t-i/Df 



(1) 



(2) 



where we restricted ourselves to the physically relevant 
situation D:$>T. 

3. From the diagrammatic language one can derive a 
formally exact kinetic equation for p{t) 



p[t) + iLop{t) 



dt' E(i - t')p{t') , 



(3) 



where Yi{t — t') is a superoperator acting on p{t') and 
is defined by the sum over all irreducible diagrams. In 
Laplace space we get the formal solution p{z) = n(z)po 
with n(z) = i/{z — Lq — iE(z)). 

4. The kernel Tj{z) is calculated by a renormaliza- 
tion group procedure. Short time scales of ^{t) are in- 
tegrated out first by introducing a short-time cutoff t^ 
into the correlation function 7(t) — > ^{t^tc)- In each 
renormalization group step, the time scales between tc 
and tc + dtc are integrated out, starting from = 
and ending aX tc = 00. As a consequence, one gener- 
ates RG equations for Lq, Gp , and the two bound- 
ary vertex operators Ap and BP (defined as the right- 
most and leftmost vertex of the kernel S(z)). Within the 



scheme of a perturbative RG analysis, generation of mul- 
tiple vertex-operators is neglected here, as in Ref. [ p2| . 
For the present model, we choose the cutoff dependence 
either as 7(t,tc) = l{t)'d{t — tc) (choice I) or 



tc) , (4) 



7(i, tc) = j^[ R{t)Q{t -tc))+ iS{t)Q{t 



with R{t) = {d/dt)R{t) (choice II). It turns out that 
choice I is better for tc ^ A, whereas choice II is needed 
for small tc in order to avoid unphysical linear dependen- 
cies on the cutoff D. At we use a smooth crossover 
from choice I to choice II. Solving the RG equations nu- 
merically gives the kernel S(z) in Laplace space, and the 
reduced density matrix can be deduced. The static sus- 
ceptibility follows from = —id/de)TTQ<TzPst, where Tro 
is the trace over the local system, and Pst denotes the sta- 
tionary solution Pst — lirat-^oo p(t) . 

5. To calculate correlation functions of the form 
x{t) — ^Ti{[az{t), (Jz]p'^'^} wc have to generalize the RG 
procedure. With C = G+ + we get in Laplace space 



X{z) - Tro{a,n(z)(C + ^c{z))Pst} , 



(5) 



where Sc(^) is analogously defined to S(z) but contains 
the "vertex" C at any time point. Within the framework 
of Ref. [Q, the RG equation for ^c{z) can be easily 
derived and reads 



_d_ 

dtc 



/>oo />f 

^c{z) = / dt I dt' 
, Jo Jo 



dt, 



Jz{t-t') j^p^-iLo{t-t')(j^-iLot' j^p' ^ ^g-j 



where BP = BP\,=o- 

Results. Figs. ^ and || show the time evolution of p{t) 
for the unbiased and biased case. Initially, the two-level 
system is prepared in the spin- up state [u). For e = 0, 
we achieve a very good agreement with CSQD. Only the 
real parts of the nondiagonal elements, which correspond 
to (ua,), exhibit a deviation of approximately 5%. How- 
ever, the CSQD can not give an accurate error for (ctj.) 
p3| . The diagonal elements oscillate in time with the 
frequency - A(A/i:))"/(i-"). This agrees with the 
renormalized tunnel matrix element, which is the char- 
acteristic energy scale of the spin boson model 0. In 
Fig. 0we obtain A^ = 0.626A (CFT gives A^ = 0.625A). 
In contrast, the real part of the nondiagonal elements is 
purely decaying for e = (only for e 7^ it also exhibits 
oscillations). In the scaling limit, i.e. for _D,A 00 
such that Ar = const, the diagonal elements are uni- 
versal, i.e. they only depend on Aj.t. The nondiagonal 
elements however have an extra factor of A^/A (see inset 
in Fig. 1^) which is consistent with |^] . The long-time be- 
haviour is given by an exponential decay. For e = the 
decay constants of the diagonal (nondiagonal) elements 
correspond to the relaxation (dephasing) time r'"'^' {t'''^p), 
see insets in Fig. 0. For t'"'^' we find good agreement with 
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CFT.For e 7^ the diagonal (nondiagonal) elements also 
contain an incoherent (coherent) decay with t'^'^p (r'"'^'). 
We note that our zero-temperature results at finite bias 
and the dephasing time are the first ones presented in 
the literature. 

The static susceptibility for the unbiased and biased 
case is shown in Figs. |3| and ^. Clearly, the susceptibil- 
ity depends strongly on the cutofi' D (see right inset in 
Fig. ^). However, in the scaling limit the susceptibility 
is universal. We normalize our results using the zero- 
temperature value Xo{T = 0)- For e = 0, the susceptibil- 
ity is a monotonous function, whereas a local maximum 
is obtained for the biased case. A quantitative compar- 
ism with recent Bethe ansatz results for the Kondo model 
shows a very good agreement for the unbiased case (see 
left inset in Fig. ||). However, for e ^ 0, there are devia- 
tions which demonstrates that the commonly used rela- 
tion of the parameters of the Kondo model to the ones 
of the spin-boson model has to be taken with care. In 
the high-temperature limit we obtain the correct (1/2T) 
law, which is independent of a, e and D |p^ . 

The result for the imaginary part of the dynamic 
susceptibility = J^^dtu e^'^^^Tr {[a S),'^z]p'"'} is 

shown in Fig. |, with S{lj) = x"{^)/^- The Shiba- 
relation provides an exact equation for the low frequency 
limit of x" at T = H 

\imSiuj) = 27raxl. (7) 

We tested Eq. (||) for different a, e and D, see Tab. |. 
For a ^ 0.1, we achieved a very good agreement (error 
smaller than 5%). Like the static susceptibility, the cor- 
relation function depends strongly on the cutoff D, but 
in the scaling limit, the normalized quantity S{lu)/S{0) 
depends only on ui/uJmax (see inset in Fig. ||), where 
^max ~ deuotcs the frequency where 5'(w) is max- 
imum. We note that our results for the correlation func- 
tion are the first ones presented for the spin-boson model 
for a = 0.1. NRG results [Hare very accurate for low 
frequency but fail for uj ~ A [Il2[ , fiow equation methods 
have already an error of ^ 25% concerning the Shiba- 
relation |l^] , and CSQD does not provide a check of the 
Shiba-relation and only gives data for the unbiased case 
e = 0(ll). 

In summary, we investigated the spin-boson model for 
ohmic dissipation using real-time renormalization group. 
For coupling parameters a ^ 0.1 — 0.2 we achieved a 
full solution of the static and dynamical properties. We 
calculated the real-time evolution of all matrix elements 
of the reduced density matrix together with the static 
susceptibility and correlation functions. In contrast to 
previous works we are not restricted to zero bias or the 
scaling limit. The restriction in a is due to the fact that 
we neglected double vertex diagrams. Our results show 
that the RTRG method is a very fiexible tool to treat 
various kinds of dissipative quantum systems. 
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FIG. 1. Real part of p{t) for a = 0.1, e = 0, D = lOOA 

and r = 0. Solid lines: RTRG. Dashed lines: CSQD; Inset: 
Rescaled real part of p{t). Solid lines: D = lOOA. Dashed 
lines: D = lOOOA. 
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FIG. 2. Real part of p(t) for a = 0.1, e = O.IA, D = lOOA 
and T = 0; Left inset: r*"^' as a function of a for D = lOOA 
and T = 0. Solid line: e = 0. Dashed line: e = 0.5A. Long 
dashed line: CFT, e = 0; Right inset: t'^"'' as a function of 
a for D = lOOA and T = 0. SoUd Une: e = 0. Dashed Une: 
e = 0.5A. 
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FIG. 3. Static susceptibility as a function of temperature 
for e = 0, D > A, and a = 0.01,0.05,0.1,0.125,0.2 (from 
top to bottom); Left inset: a = 0.125. Solid line: RTRG. 
Dashed line: Bethe ansatz; Right inset: Cutoff dependence 
for a = 0.1. Solid line: D = lOOA. Dashed line: D = lOOOA. 
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FIG. 4. Static susceptibility as a function of temperature 
for a = 0.1, L> > A, and e/A = 0, 0.1, 0.2, 0.5, 1 (from top to 
bottom). 
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FIG. 5. 5 as a function of a; for a = 0.1, T = 0. 

Solid line: e = 0, D = lOOA. Dotted line: e = O.IA, 
D = lOOA. Dashed line: e = 0, D = lOOOA. Dot-dashed 
line: e = O.IA, D = lOOOA; Inset: Rescaled S. Upper curve: 
e = 0. Lower curve: e = O.IA (for D = lOOOA we rescaled 
e = (1000/100)"/^^-"' X O.IA = 10^/^ x O.IA). 



a 


e/A 


D/A 


XoA 


lim A^S(uj) 


error 


0.01 


0.0 


100 


1.0511 


0.0680 


2.04% 


0.05 


0.0 


100 


1.2899 


0.5083 


2.80% 


0.1 


0.0 


100 


1.6868 


1.8343 


2.57% 


0.2 


0.0 


100 


3.2400 


12.3085 


6.93% 


0.01 


0.0 


1000 


1.0772 


0.0711 


2.50% 


0.05 


0.0 


1000 


1.4534 


0.6445 


2.92% 


0.1 


0.0 


1000 


2.1804 


3.0880 


3.33% 


0.2 


0.0 


1000 


5.7219 


35.4919 


14.75% 


0.1 


0.01 


100 


1.6859 


1.8309 


2.49% 


0.1 


0.05 


100 


1.6671 


1.7523 


0.34% 


0.1 


0.1 


100 


1.6036 


1.6209 


0.32% 


0.1 


0.01 


1000 


2.1778 


3.0911 


3.66% 


0.1 


0.05 


1000 


2.1369 


2.7233 


5.21% 


0.1 


0.1 


1000 


2.0187 


2.4377 


4.91% 



TABLE I. Shiba-relation for different a, e and D. 
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